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Macroscopic material properties from quasi-static, microscopic 
simulations of a two-dimensional shear-cell 

Marc Latzel, Stefan Luding, and Hans J. Herrmann (*) 



Abstract One of the essential questions in the area 
of granular matter is, how to obtain macroscopic tenso- 
rial quantities like stress and strain from "microscopic" 
quantities like the contact forces in a granular assembly. 
Different averaging strategies are introduced, tested, and 
used to obtain volume fractions, coordination numbers, 
and fabric properties. We derive anew the non-trivial re- 
lation for the stress tensor that allows a straightforward 
calculation of the mean stress from discrete element simu- 
lations and comment on the applicability. Furthermore, we 
derive the "elastic" (reversible) mean displacement gradi- 
. ent, based on a best-fit hypothesis. Finally, different com- 
binations of the tensorial quantities are used to compute 
some material properties. 

The bulk modulus, i.e. the stiffness of the granulate, 
is a linear function of the trace of the fabric tensor which 
itself is proportional to the density and the coordination 
number. The fabric, the stress and the strain tensors are 
not co-linear so that a more refined analysis than a clas- 
sical elasticity theory is required. 

; 1 

Introduction 

Macroscopic continuum equations for the description of 
the behavior of granular media rely on constitutive equa- 
tions for stress, strain, and other physical quantities de- 
scribing the state of the system. One possible way of ob- 
taining an observable like the stress is to perform discrete 
particle simulations Jl|, |) and to average over the "micro- 
scopic" quantities in the simulation, in order to obtain an 
averaged macroscopic quantity. In the literature, slightly 
different definitions for stress and strain averaging proce- 
dures can be found §-(l§. 
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The outline of this study is as follows. In section g 
the discrete element simulation method is discussed and, 
in section [| some averaging methods are introduced and 
applied to a scalar quantity, namely the volume fraction. 
Section ^ contains the definitions and averaging strategies 
for fabric, stress, and elastic strain and in section ^ some 
material properties are extracted from the results obtained 
in section [|. 

2 

Modelling discrete particles 

The elementary units of granular materials are mesoscopic 
grains. In order to account for the excluded volume, one 
can assume that the grains are impenetrable but deform 
under stress. Since the realistic modelling of the defor- 
mations of the particles in the framework of a continuum 
theory Jl^,|lJ] would be much too complicated, we relate 
the interaction force to the overlap 8 of two particles. Note 
that the evaluation of the inter particle forces based on the 
overlap may not be sufficient to account for the nonlinear 
stress distribution inside the particles. Consequently, our 
results presented below are of the same quality as this 
basic assumption. 

The force laws used are material dependent, involving 
properties such as Young's modulus of elasticity, and have 
to be validated by comparison with experimental measure- 
ments |J|-|l7|]. 

When all forces /; acting on the particle i, either from 
other particles, from boundaries or from external forces, 
are known, the problem is reduced to the integration of 
Newton's equations of motion for the translational and 
the rotational degrees of freedom 

d 2 d 2 
mi— ri = , and h-^&i = Mi . (1) 

The mass of particle i with diameter di is m,, and its 
moment of inertia is Ii = qimi(ai) 2 , with the radius = 
di/2 and the dimensionless shape factor qi. The vectors 
r*i and give the position and the orientation angles 
of particle i, respectively. In our model attractive forces 
and the presence of other phases are neglected, we focus 
on "dry granular media" . Particle-particle interactions are 
short range and active on contact only. The total force 
(torque) due to other particles is thus f t = fl (Mi = 
J2 C M^), where the sum runs over all contacts of particle 
i. The torque M.\ = if X is related to the force f\ 
via the cross product with the branch vector Z£ from the 
particle center to the point of contact c. Eq. (|l|) consists 
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of six scalar equations in three dimensions and reduces 
to three equations in two dimensions (2D), two for the 
linear and one for the rotational degree of freedom. In the 
following the force laws for f\ accounting for excluded 
volume, dissipation, and friction are introduced. 



2.1 

Force laws 

The particles i and j interact only when they are in con- 

-TV, I 



tact so that their overlap S = \((k + dj) — (r 
is positive, with the unit vector n 
that points from j to i. The symbol '•' denotes the scalar 
product of vectors or, more generally, the order-reduction 
by one for each of two tensors. 

The first contribution to the force acting on particle i 
from j is an elastic repulsive force 



/n,cl = k n 5n 



(2) 



where k n is proportional to the material's modulus of elas- 
ticity with units [N/m]. Since we are interested in disks 
rather than spheres, we use a linear spring that follows 
Hooke's law, whereas in the case of elastic spheres, the 
Hertz contact law would be more appropriate jl8], |l^] . 

The second contribution, a viscous dissipation, is given 
by the damping force in the normal direction 



/n,diss = InSn , 



(3) 



where 7„ is a phenomenological normal viscous dissipation 
coefficient with units [kg s^ 1 ] and S — —Vij ■ n the relative 
velocity in the normal direction w, 3 = Vi — Vj . 

The third contribution to the contact force - account- 
ing for tangential friction - can be chosen in the sim- 
plest case, according to Coulomb, as /t, friction < — A*|/n|*j 
where /i is the friction coefficient and t = £/|£| is the tan- 
gential unit- vector parallel to the tangential component of 
the relative velocity £ = Vij — (v^ -n)n. Because this non- 
smooth ansatz leads to numerical problems for small £ a 
regularizing viscous force /t, viscous = — 7t£ is added. The 
two forces are combined by taking the minimum value 



ft = -min(|7ii|, \fif n \)t 



(4) 



The effect of a more realistic tangential force law according 
to Cundall and Strack Q wm be reported elsewhere ||2l| . 
Due to the boundary conditions introduced below, it is 
also necessary to account for the friction with the bottom 



fb 



-fibmgv 



(5) 



with the unit vector in the direction of the particles veloc- 
ity v = vj\v\ . The effect of a bottom friction was discussed 
also in B2l . In summary, we combine the forces and use 



fi 



(,fn,e\ H~~ ^n,diss + ft) + f b 



(6) 



for the forces acting on particle i at its contact c with 
particle j. 



2.2 

Model system 

In the simulations presented in this study, a two-dimensio- 
nal Couette shear-cell is used, filled with a bidisperse pack- 
ing of disks, as sketched in Fig. [l]. The system undergoes 
slow shearing introduced by turning the inner ring. The 
properties of the particles, used for the force laws of sub- 
1 , are summarized in table [lj The boundary con- 



section 



2.1 



a) filling 




Fig. 1. A schematic plot of the model system 



ditions are based on an experimental realization [|23| - 
For more details on other simulations, see 122, E4L 



property 


values 


diameter d sma ,n, mass m sma n 


7.42 mm, 0.275 g 


diameter di ar g C , mass mi arge 


8.99 mm, 0.490 g 


wall-particle diameter d wa n, 


2.50 mm 


system/disk-height h 


6 mm 


normal spring constant fc n 


352.1 N/m 


normal viscous coefficient 7„ 


0.19kg/s 


tangential viscous damping 7t 


0.15kg/s 


Coulomb friction coefficient fi 


0.44 


bottom friction coefficient /ib 


2 x 10" 5 


material density po 


1060 kgm" 3 



Table 1. Microscopic material parameters of the model. 



In the simulations different global volume fractions 



1 



N 

E 

P =i 



V" 



(7) 



of the shear-cell are examined. The sum in Eq. (Q) runs 
over all particles p with volume V p in the cell with I4ot = 
tt(RI — In this study v is varied between 0.8084 and 
0.8194, by varying the particle number, see table |[ For 
the calculation of the global volume fraction, the small 
particles glued to the wall are counted with half their vol- 
ume only, and thus contribute with r? wa n = 0.0047 to v. 

Note that we use the phrases "volume" and "volume 
fraction" even if, strictly speaking, the unfamiliar terms 
"disk area" and "area fraction" could be used. The reasons 
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global volume 


number of particles 


^max 


At 




fraction v 


small 


large 


(s) 


(») 


A 


0.8084 


2511 


400 


335 


5 


B 


0.8103 


2520 


399 


119 


1 


C 


0.8133 


2524 


404 


119 


1 


D 


0.8149 


2545 


394 


119 


1 


E 


0.8180 


2538 


407 


505 


5 


F 


0.8194 


2555 


399 


119 


1 



Table 2. Details of the simulation runs studied in this paper. 



for this choice are: (i) The methods discussed in this study 
are straightforwardly generalized to three dimensions and 
(ii) the particles are three dimensional objects with height 
h anyway, so that the use of the word "volume" is justified. 



2.3 

Initial conditions and steady state 

The simulations are started in a dilute state with an ex- 
tended outer ring R a (t = 0) > R — 0.2524 m, and the 
inner ring already rotates counterclockwise with constant 
angular frequency Q = 2-njTi = 0.1 s -1 and period 7j = 
62.83 s. The extended outer ring is used in order to allow 
for a random, dilute initial configuration. The desired den- 
sity is then reached by reducing the volume. The radius 
of the outer ring is reduced within about two seconds to 
reach its desired value R . Afterwards, the outer ring is 
kept fixed and the inner ring continues to rotate until at 
t = tmax the simulation ends. If not explicitly mentioned, 
averages are performed after about one rotation at t — 60 s 
(to get rid of the arbitrary initial configuration) , and dur- 
ing about one rotation, until t = 119 s. 



From the micro- to a macro-description 

In the previous section, the microscopic point of view was 
introduced, as used for the discrete element method. Par- 
ticles are viewed as independent entities which interact 
when they come in contact. In this framework, the know- 
ledge of the forces at each contact is sufficient to model the 
dynamics and the statics of the system. Tensorial quan- 
tities like the stress or the deformation gradient are not 
necessary for a discrete modelling. However, subject of 
current research is to establish a correspondence to con- 
tinuum theories by computing tensorial quantities like the 
stress <r, the strain e, as well as scalar material properties 
like, e.g., the bulk and shear moduli [0,^,[l^]. In the course 
of this process, we first discuss averaging strategies using 
the material density as an example. 



3.1 

Averaging strategies 

Most of the measurable quantities in granular materials 
vary strongly on short distances. Thus, computing aver- 
ages necessitates dealing with or smearing out the fluctu- 
ations. In order to suppress the fluctuations, we perform 
averages in both time and space. This is possible due to 



the chosen boundary condition. The system can run for 
long time in a quasi-steady state and, due to the cylin- 
drical symmetry, points at a certain distance R from the 
origin are equivalent to each other. Therefore, averages 
are taken over many snapshots in time with time steps At 
and on rings of material at a distance r = R — Ri from 
the inner ring. The width of the averaging rings is Ar, 
so that the averaging volume of one ring is V r — 2irrAr. 
For the sake of simplicity (and since the procedure is not 
restricted to cylindrical symmetry), the averaging volume 
is denoted by V = V r in the following. The rings are num- 
bered from s = to B — 1, with B = (R Q — Ri)/Ar, and 
ring s reaches from r s = r — Ar/2 to r s +i = r + Ar/2. The 
averaging over many snapshots is somehow equivalent to 
an ensemble average. However, we remark that different 
snapshots are not neces sari ly independent of each other as 
discussed in subsection 3.4. Also the duration of the simu- 



lation maybe not long enough to explore a representative 
part of the phase space. 

The cylindrical symmetry is accounted for by a rota- 
tion of all directed quantities like vectors, depending on 
the cartesian position r, = (xi,yi) of the corresponding 
particle i. The orientation of particle i is 4>i = BXctBXi{yi/xi) 
for Xi > and periodically continued for %i < so that 4>i 
can be found in the interval [— tt, tt]. The vector n c that 
corresponds to contact c of particle i is then rotated about 
the angle — 4>i from its cartesian orientation before being 
used for an average. Note that this does not correspond to 
a transformation into orthonormal cylindrical coordinates. 

Finally, we should remark that the most drastic as- 
sumption used for our averaging procedure is the fact, that 
all quantities are smeared out over one particle. Since it 
cannot be the goal to solve for the stress field inside one 
particle, we assume that a measured quantity is constant 
inside the particle. This is almost true for the density, but 
not e.g. for the stress. However, since we average over all 
positions with similar distance from the origin, i.e. aver- 
ages are performed over particles with different positions 
inside a ring, details of the position dependency inside the 
particles will be smeared out anyway. An alternative ap- 
proach was recently proposed by I. Goldhirsch [jll| who 
smeared out the averaging quantities along the lines con- 
necting the centers of the particles and weighed the con- 
tribution according to the fraction of this line within the 
averaging volume. 



3.2 

Volume fractions 

The first quantity to measure is the local volume fraction 



(8) 



with the particle volume V p . w v is the weight of the par- 
ticle's contribution to the average. 

This formalism can be extended to obtain the mean 
value of a quantity Q in the following way: 



(9) 



p&V 



4 



with the pre-averaged particle quantity 

c p 

q p = J2Q c - (io) 

c-1 

with the quantity Q c attributed to contact c of particle p. 
In the following the brackets (Q) denote the average over 
a quantity Q of the particles in an averaging volume V. 

The simplest choice for w v is to use w v — 1, if the 
center of the particle lies inside the ring, and w v = oth- 
erwise. This method will be referred to as particle- center 
averaging in the following. 

A more complicated way to account for those parts of 
the particles which lie partially inside a ring is to use only 
the fraction of the particle volume that is covered by the 
averaging volume. Since an exact calculation of the area of 
a disk that lies in an arbitrary ring is rather complicated, 
we assume that the boundaries of V are locally straight, 
i.e. we cut the particle in slices, as shown in Fig. ||. This 
method will be referred to as slicing in the following. The 
error introduced by using straight cuts is well below one 
per-cent in all situations considered here. 




Fig. 2. Schematic plot of a particle p at radial position r v 
which is cut into pieces by the boundaries r s of the averaging 
volumes. We assume s = 0, . . ., m + 1 such that all r s with 
s = 1, . . ., m hit the particle, i.e. \r v — r s \ < d/2. 

The volume Vf = WyV p of a particle p which partially 
lies between r s and r s+ i is calculated by subtracting the 
external volumes V~ and from the particle volume 
V p = n{d/2) 2 so that 

v* = vp- v- - v+ 

= {d/2) 2 [n -<t> s + sin(> s ) cos(0 s ) 

- s+ i + sin(0 s+ i) cos(0 s+ i)] (11) 

with <j) s = arccos(2(r p — r s )/d) and <j) s +i = arccos(2(r s+ i — 
r p )/d). The term (d/2) 2 <p is the area of the segment of the 
circle with angle 2cf>, and the term (d/2) 2 sin(</>) cos(0) is 
the area of the triangle belonging to the segment. In Fig. || 
the case s = 1 is highlighted, and the boundaries between 
V s ~~ , V[, and V, + are indicated as thick solid lines. The 
two outermost slices Vq = V{~ and VJ = 1 have to be 
calculated separately. 

When the two averaging strategies (particle-center and 
slicing) are compared for different widths Ar of the aver- 
aging rings, so that the number of intervals is B — 20, 40, 
or 60. For Ar « d S maii (B = 20), the results are almost in- 
dependent of the averaging procedure. For finer binning 



B s > 30, the slicing procedure converges on a master 
curve with weak oscillations close to the walls which are 
due to the wall-induced layering of the particles. The re- 
sults of the particle-center averaging strongly fluctuate for 
B c > 24. These oscillations come from ordered layers of 
the particles close to the walls so that the slicing method 
reflects the real density distribution for fine enough bin- 
ning. The particle-center method, on the other hand, leads 
to peaks, where the centers of the particles in a layer are 
situated and to much smaller densities where few particle 
centers are found; the particle-center density is obtained 
rather than the material density. 



3.3 

Representative elementary volume REV 

An important question is, how does the result of an av- 
eraging procedure depend on the size of the averaging 
volume V. We combine time- and space averaging, i.e. 
we average over many snapshots and over rings of width 
Ar, so that the remaining "size" of the averaging vol- 
ume is the width of the rings Ar. In Fig. || data for v 
at fixed position R = 0.12, 0.13, 0.14, and 0.20m, but 
obtained with different width Ar, are presented. The po- 
sitions correspond to r/d sma n « 2.2, 3.6, 4.9, and 13, when 
made dimensionless with the diameter of the small parti- 
cles. Both the particle-center method (open symbols) and 
the slicing method (solid symbols) are almost identical 
for Ar/d STaa \\ > 2, for the larger Ar the averaging vol- 
ume can partially lie outside of the system. For very small 
Ar/d svaa \\ < 0.1 the different methods lead to strongly dif- 
fering results, however, the values in the limit Ar — ► are 
consistent, i.e. independent of Ar besides statistical fluc- 
tuations. In the intermediate regime 0.1 < Ar/d sma \\ < 2, 
the particle-center method strongly varies, while the slic- 
ing method shows a comparatively smooth variation. 

Interestingly, all methods seem to collapse at Ar^EM 
~ ^smaii (and twice this value), nearly to the size of the 
majority of the particles. For the examined situations, we 
observe that the particle-center and the slicing method 
lead to similar results for 0.97 < <4rREv/^smaii < 1-03. 
This indicates that the systems (and measurements of 
system quantities) are sensitive to a typical length scale, 
which is here somewhat smaller than the mean particle 
size. When using this special Ar^y value, one has B = 20 
or B = 21 binning intervals. The open question of this be- 
ing a typical length scale that also occurs in systems with 
a broader size spectrum, cannot be answered with our 
setup, due to the given particle-size ratio. 



3.4 

Time averaging 

In order to understand the fluctuations in the system over 
time, and to test whether subsequent snapshots can be 
assumed to be independent, the volume fraction v is dis- 
played for snapshots at different averaging times in Fig. 

Changes in density are very weak and mostly occur in 
the dilated shear zone for small r. From one snapshot to 
the next, we frequently find, that the configuration in the 
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Fig. 3. Volume fraction v at different distances r from the 
inner ring, plotted against the width Ar of the averaging ring 
from simulation D, see table ^. Note that the horizontal axis 
is logarithmic. The open symbols are results obtained with the 
particle-center method, the solid symbols are results from the 
slicing method. The inset is a zoom into the large Ar region, 
and the arrows indicate the optimal width Ar for the particle- 
center method for which the results are almost independent of 
the averaging procedure. 
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Fig. 4. Volume fraction u, plotted against the dimensionless 
distance r/d ama ii from the inner ring, for different times from 
simulation E. Open symbols belong to six subsequent snap- 
shots with At — 5 s, the small, solid symbols are snapshots 
after longer times. 



outer part of the shear cell has not changed, whereas a 
new configuration is found in the inner part. Only after 
rather long times does the density change also in the outer 
part. Thus, simulation results in the outer part are sub- 
ject to stronger fluctuations because the average is per- 
formed over less independent configurations than in the 
inner part. 



Macroscopic tensorial quantities 

In this section, the averaged, macroscopic tensorial quan- 
tities in our model system are presented. The fabric tensor 
describes the contact network, the stress tensor describes, 
in this study, the stress due to normal forces, and the 
deformation gradient is a measure for the corresponding 
elastic, reversible deformations. A more detailed analysis, 
where the tangential forces are also taken into account, is 
in progress plj ] and some details in that direction can be 
found in Refs. |,[To[|l2|. However, we checked the influ- 
ence of the tangential forces in our system and found their 
effect to be always smaller than ten per-cent. 



4.1 

Micro-mechanical fabric tensor 

In assemblies of grains, the forces are transmitted from 
one particle to the next only at the contacts of the parti- 
cles. In the general case of non-spherical particles, a pack- 
ing network is characterized by the vectors connecting the 
centers of the particles and by the particle-contact vec- 
tors. Furthermore, the local geometry of each contact is 
,27 1, see Fig. ||. In our case, with spherical 



important 

particles, the situation is simpler with respect to both the 
spherical contact geometry and the fabric. 




Fig. 5. Schematic plot of a particle with radius a and C p = 4 
contacts as indicated by the small circles. The branch vector 
l pc and the normal unit vector n c are displayed at contact 
c = 1. 



4.1.1 

Fabric tensor for one particle 

One quantity that describes the contact configuration of 
one particle to some extent, is the second order fabric 
tensor \M,WH 



c=l 



(12) 



where n c is the unit normal vector at contact c of parti- 
cle p with C p contacts. The symbol ® denotes the dyadic 
product in this study. Other definitions of the fabric use 
the so-called branch vector l pc from the center of parti- 
cle p to its contact c, however, the unit normal and the 
unit branch vector are related by a p n = l pc in the case of 
spheres or disks, so that the definition 



1 



c p 



,2 Z-^i 



l pc ® l pc 



(13) 
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Fig. 6. a-c Plot of (a) Fv — tr(F), (b) Fd/Fv = dev(F)/tr(F), and (c) 4>f against the dimensionless distance r/d sma jx from 
the inner ring. The global volume fraction is given in (b) and valid also for (a) and (c). 



is identical to Eq. (p"2|). The fabric tensor in Eq. fll2| ) is 
symmetric by definition and thus consists of up to three 
independent scalar quantities in two dimensions. The first 
of them, the trace (or volumetric part) Fy — tr(F p ) = 
{Fmax + -Fmin), is the number of contacts of particle p, 
with the major and the minor eigenvalues -Fmax and F m { n , 
respectively. One gets from Eqs. (|lj) or ( |l3| ) the number 
of contacts of particle p 

c p 

tr(F p ) = tr ( nC ® nC ) = CP > ( 14 ) 

c=l 

since the scalar product of n c with itself is unity by defi- 
nition. The second scalar, the deviator Fd — F max — F m ; n , 
accounts for the anisotropy of the contact network to first 
order, and the third, the angle <pF, gives the orientation of 
the "major eigenvector" , i.e. the eigenvector correspond- 
ing to F max , with respect to the radial direction. In other 
words, the contact probability distribution is proportional 
to the function [F v + F D cos{2(<p - </> F ))] @||,[29), when 
averaged over many particles, an approximation which is 
not always reasonable ]3C[ |. 



4.1.2 

Averaged fabric tensor 

Assuming that all particles in V contribute to the fabric 
with a weight of their volume V p one has 

c p 

pgV c =l 

The values of Fy, Fp/Fy, and <f>p, as obtained from our 
simulations, are plotted in Figs. ^(a-c). The trace of the 
fabric tensor and thus the mean number of contacts in- 
creases with increasing distance from the inner ring, and is 
reduced in the vicinity of the walls. With increasing mean 
density, the trace of F is systematically increasing, while 
the deviatoric fraction seems to decrease; this means that 
a denser system is slightly more isotropic concerning the 
fabric. The major eigendirection is tilted counterclockwise 
by somewhat more than 7r/4 from the radial outward di- 
rection, except for the innermost layer and for the strongly 
fluctuating outer region. 



In analogy to the trace of the fabric for a single parti- 
cle, the trace of the averaged fabric is 

tr(F) = i^<F^ = (^) , (16) 
pev 

which, in the case of a regular, periodic contact network of 
almost identical particles with a p ~ a, reduces to the sum 
over all of particles in V with the prefactor v defined in 
Eq. (||). Now, one has a relation between the coordination 
number, i.e. the mean number of contacts per particle C = 
C(r) = {C p )/v : the volume fraction v, and the averaged 
fabric F, as a combination of Eqs. (|J) and (|l6|): 

tr(F) ~ vC . (17) 

As a test for the averaging procedure, we plot in Fig. [7] 
tr(F) against vC and obtain all data points from all sim- 
ulations collapsing close to the identity curve. For a theo- 
retical derivation of the small (about 1 per-cent) deviation 
due to the polydispersity, see Ref. Q. 
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4.2 

Micro-mechanical stress tensor 

The micro-mechanical stress tensor is derived in a way 
similar to Ref . || . For an arbitrary volume V with surface 
dV, the mean stress is defined as 



(18) 



where er' = er'(x) is the position dependent local stress 
inside V and x is the coordinate integrated over. Note 
that er' can be a strongly fluctuating function of x. 

If the stress in the pore-space vanishes, the integral 
can be split into a sum over the stresses er p , pre-averaged 
for particles p, with the respective center of mass vectors 
x p . The mean stress is thus 

CT = < CT > = v £ w v VP(tP = v £ < / , dy ' <T ' ' ( 1Q ) 



which corresponds to a smearing out of er' over the parti- 
cles, with er p , the average stress in particle p, to be derived 
in the following. 

One can rewrite the transposed stress tensor so that 



er T = grad x ■ cr T = div [x ® er) — x ® div er 



(20) 



by introducing the unit tensor I = grad x, and by applying 
the series rule in the first term on the right hand side. Note 
that the tensors er on the right hand side are transposed 
with respect to the left hand side er T . The transposed 
stress is used for the following operations for the sake of 
simplicity, however, using the stress directly should lead 
to the same result. For a more detailed treatment, see 
Ref. . In static equilibrium and in the absence of body 
forces, the term div er on the right hand side vanishes and 
er is symmetric. 



4.2.1 

Mean stress for one particle 

Using the definition of the transposed stress tensor in 
Eq. (20) the integral over one particle in Eq. 

1 { dV div (x (g) er') 
Jvp 



91) becomes 



VP 
1 

V p ,i dV p 



= — dS (x <g> er') • n 



(21) 



by application of Gauss' theorem. In Eq. (|T|) dS is the 
surface element of V p on dV p and n is the outwards nor- 
mal unit vector. Using the definition of the stress vector 
t = er' ■ n one arrives at 

W* = lkl dSx®t. (22) 

A force f c acting at a contact c with area Ss c leads to 
a stress vector t c = f c /Ss c - even in the limit of small 
contact area 5s c <C a p . Here, we apply the simplifying as- 
sumption that the force f c is constant on the surface 6s c , 
i.e. we do not resolve any details at the contact. Explicitly 
writing the surface element dS as a sum, 

c 

dS = ^2S(\x-x c \) 6s c , (23) 



leads, after integration over the delta functions, to the 
transposed stress 



C" 



C— 1 



(24) 



Transposing Eq. ( |24| ) leads to an exchange of x c and / c , 
and thus to the expression for the mean stress inside par- 
ticle p: 



c p 



yp 



(25) 



c=l 



which is, at first glance, dependent on the frame of refer- 
ence associated with the vector x c . 

Using vector addition, one has x c = x p + l pc , with the 
position vector x p of particle p and the branch vector l pc , 
pointing from the center of mass of particle p to contact 
c. Inserting this relation in Eq. ( ^5|) leads to 



\C— 1 / C— 1 



■pc 



VP 



(26) 



since the first sum vanishes in static equilibrium, where 

c p 



4.2.2 

Averaged stress tensor 

Inserting Eq. (fH in Eq. (|l|) gives a double sum over all 
particles with center inside the averaging volume V and 
all their contacts 



i pc 



(28) 



c=l 



In other words, stress is pre-averaged over all particles and 
then averaged over V. 

Inserting Eq. (|B|) in Eq. ( |l9|) leads to the mathemat- 
ically identical sum 



c 

P f c 



n =vl^ w 'vL^ 

pev c=i 



(29) 



This expression can be transformed into a sum over all 
contacts c(p, q) € V, where at least one of the participat- 
ing particles p and q lies inside the volume V. The force 
f c = f pq = —f qp , acting at contact c, from particle q on 
p is equal, but opposite in direction to the force acting at 
the same contact from particle p on q, so that 



q f qp\ 



(30) 
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since Wy — w v = when both particles are completely in- 
side V. Note that every contact is visited once in Eq. ( [j0|) 
but twice in Eqs. (28) and (|2^). The expression c 6 dV 
means for the particle-center averaging method that con- 
tacts contribute only if one particle is inside V while the 
other is outside. In the framework of the slicing method, 
contacts contribute if at least one particle is cut by dV so 
that w v — w v 7^ 0. 

Given an arbitrary averaging volume V, and all the 
information about the forces at all contacts from the dis- 
crete element simulations, it is obvious that - from the 
technical point of view - the sum in Eq. (28) can be done 
relatively easily, whereas the sum in Eq. (3C) requires the 
identification of the contacts at the "surface" of V, before 
the possibly shorter summation can be performed. Note 
that both expressions Eq. ( |28| ) and Eq. (f30|) have been 
used by different authors, see e.g. Jo|, |s|, plf . 

The values of the volumetric stress ay = tr(er), the 
deviatoric fraction (Jd/o~v = dev(er)/tr(cr), as obtained 
from our simulations, are plotted in Figs. ||(a-b). The vol- 
umetric stress is constant, besides fluctuations, whereas 
the deviatoric fraction is rather large in the shear-zone and 
decays with increasing distance from the inner ring, like 
the fabric. Also, the deviatoric fraction of the stress ap- 
pears slightly reduced in magnitude for increasing global 
density. Note that the stress varies over almost two or- 
ders of magnitude, while the mean density is changed only 
weakly from v = 0.8084 to 0.8194, see table |. In Fig. |= 4> a 
describes the angle which the principal axis of the stress 
tensor is tilted from the outward direction. The stress ten- 
sor is, on average, rotated counter clockwise by somewhat 
less than 7r/4 from the outward direction. It is only in the 
outermost part that strong fluctuations exist around the 
mean. We now provide evidence that the fabric and stress 
are not co-linear. 



4.3 

Mean strain tensor 

To achieve the material properties of a granular ensem- 
ble one is interested in the stress-strain relationship of 
the material. One of the simplest techniques used, is the 
application of "Voigt's hypothesis" which assumes that 
the strain is uniform and that every particle displacement 
conforms to the mean displacement field ]32]| . Thus, the 
expected displacement at contact c of particle p, relative 
to the force free situation, and due to the mean displace- 
ment gradient e, is e • l pc , with l pc = a p n c , and the mean 
contact deformation S. Note that the linear, symmetric 
strain e = |(e + e T ) is not identical to the displacement 
gradient, in general. In our study, we follow the approach 
of Liao et al. poL who used l pg instead of l pc , and as- 
sume that the actual displacement field does not coincide 
with, but fluctuates about the mean displacement field. 
The difference between the actual displacement A pc and 
the expected displacement is 



x pc = e • l pc - A pc 



(31) 



The actual displacement is directly related to the simula- 
tions via A pc — 5 c n c . 



If one assumes that the mean displacement field best 
approximates the actual displacement, one can apply a 
"least square fit" to the total fluctuations 

c 

5 =-j7E<I> PC ) 2 ' ( 32 ) 
P ev c=i 

by minimizing S with respect to the mean displacement 
gradient so that 



dS_ 







V ^ 



(33) 



C p 

<E 

c=l 



(e • l pc 



A pc ) ■ —(e-l pc 
oe 



A pc ) 



These four equations for the four components of e in 2D 
can be transformed into a relation for the mean displace- 
ment tensor as a function of the contact displacements and 
the branch vectors. By assuming that dA pc /de — 0, the 
expression "-^(e-Z pc )" becomes a dyadic product "®Z pc ", 
as can be seen by writing down the equation in index no- 
tation. Extracting e from the sum (what leaves l pc (g> l pc , 
the core of the fabric, in the first term) and multiplying 
the equation with A = F _1 , the inverse of the fabric, see 
Eq. (pi), we find that 



pgV c=l 



• l pc 



A . 



(34) 



The values of ey, en/ey, and <fi e , as obtained from our 
simulations, are plotted in Figs. |^(a-c). The elastic, volu- 
metric deformation gradient of the granulate is localized 
in the shear zone and the effect is stronger in the less dense 
systems. Due to dilation it is easier to compress the dilute 
material closer to the inner ring compared to the outer 
part. Like the fabric and stress, the strain also becomes 
more isotropic with increasing mean density. 



5 

Results 

At first, we compute the mean-field expectation values for 
<j and e, in order to get a rough estimate for the orders of 
magnitude of the following results. Replacing, in Eq. (p8|), 
f° by its mean / = k n 5, a p by a = a, and l pc by cm c , one 
gets 



er = (k„ S/na) F 



(35) 



Performing some similar replacements in Eq. (p4j) , leads 
to 

e = (8 1 a) I , (36) 



e = (n/k n ) tr • A . (37) 

The material stiffness, E, can be defined as the ratio of the 
volumetric parts of stress and strain, so that one obtains 
from Eq. @ and @ 



E = (fc„/27r) tr(F) 



(38) 



In Fig. |lCj the rescaled stiffness of the granulate is plotted 
against the trace of the fabric for all simulations. Note 
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that all data collapse almost on a line, but the mean-field 
value underestimates the simulation data by a few per- 
cent. Simulation data for different k n and even data from 
simulations with neither bottom- nor tangential friction 
collapse with the data for fixed k n and different volume 
fractions, shown here. The deviatoric fraction of <x is plot- 
ted in Fig. [ll] against the deviatoric fraction of F. Here, 
the data are strongly underestimated by the mean-field 
result in Eq. (|35|). A similar plot for the deviatoric frac- 
tion of the deformation gradient against Fd / Fy shows a 



gathering of the data close to the identity curve, in dis- 
agreement with Eq. (pf^). Thus, we conclude that the devi- 
atoric parts of stress, deformation gradient and fabric are 
interconnected in a more complicated way than suggested 
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by the simple mean field estimates [Ql |32|, 

In Fig. 1^ the ratio of the deviatoric parts of stress and 
strain is plotted against the trace of the fabric. Like the 
material stiffness, both quantities are proportional, only 
G shows much stronger fluctuations and has a proportion- 
ality factor of about 1/3. We did not use the traditional 
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Fig. 12. Scaled granulate shear resistance G = dev(er)/dev(e) 
plotted against tr(F) from all simulations. 

definition of the shear modulus since our tensors are 
not co-linear as shown below. 

In Fig. [l^ the orientations of the tensors are plotted 
for some simulations and in the inner part of the shear 
cell. In the outer part, the deviatoric fraction is usually 
around 10 per-cent, i.e. so small that the orientations be- 
come noisier. Simulations A and B are skipped here, be- 
cause the data of the orientations are rather noisy due to 
the low density which leads to intermittent behavior with 
strong fluctuations. We observe that all orientation angles 
cj) show the same qualitative behavior, however, the fabric 
is tilted more than the stress which in turn is tilted more 
than the deformation gradient, where the orientations are 
measured in counter clockwise direction from the radial 
outward direction. 

As a final cross-check, inserting the measured values 
for a and A into Eq. ( |37| ) leads to the measured values of 
e within a few percent deviation. However, the orientation 
of the deformation gradient is not well reproduced - it 
seems very sensitive to small fluctuations. 
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Fig. 13. Orientation of the tensors F, a, and e plotted against 
r/dsmaii from simulations C, D, E, and F with B s — 60. The 
data are shown in the range < r/d sma n < 10 only. Open 
symbols are fabric, open symbols connected by lines are stress, 
and solid symbols are strain, the mean densities corresponding 
to the symbols is given in the plot. 
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Summary and Conclusion 

Discrete element simulations of a 2D Couette shear cell 
were presented and used as the basis for a micro-macro 
averaging procedure. In the shear cell a shear band is lo- 
calized close to the inner, rotating cylindrical wall. The 
boundary conditions were chosen to allow for averaging 
over large volumes (rings with width Ar) and over a steady 
state and thus over long times. The configurations changed 
rather rapidly in the shear band, whereas the system is 
frozen in the outer part, a fact which requires either ex- 
tremely long simulations or a sampling over different ini- 
tial conditions in order to allow an ensemble averaging in 
the outer part. 

The simplest averaging strategy involves only the par- 
ticle-centers as carriers of the quantities to be averaged 
over, whereas a more advanced method assumes the quan- 
tities to be homogeneously smeared out over the whole 
particle which is cut in slices by the averaging volumes. 
Both methods agree if the averaging volumes are of the 
particle size (or multiples), but for other sizes differing 
results are obtained. The slicing method shows discretiza- 
tion effects in the range of averaging volume widths from 
one to one fifth of a particle diameter, while the particle- 
center method shows fluctuations due to the choice of the 
averaging volume in a much wider range. 

The material density, i.e. the volume fraction, the co- 
ordination number, the fabric tensor, the stress tensor and 
the elastic, reversible deformation gradient were obtained 
by the averaging procedures. The fabric is linearly propor- 
tional to the product of volume fraction and coordination 
number. In the shear band, dilation together with a re- 
duction of the number of contacts is observed. The mean 
stress is constant in radial direction while the deformation 
gradient decays with the distance from the inner wall. The 
ratio of the volumetric parts of stress and strain gives the 
effective stiffness of the granulate, which is small in the 
shear band and larger outside, due to dilation. 

In the shear band, large deviatoric components of all 
tensorial quantities are found, however, decreasing with 
increasing distance from the inner wall. The isotropy of the 
tensors grows only slightly with increasing density and all 
tensors are tilted counterclockwise from the radial direc- 
tion by an angle of the order of 7r/4. The system organizes 
itself such that more contacts are created to act against 
the shear. An essential observation is that the macroscopic 
tensors are not co-linear, i.e. their orientations are differ- 
ent. The orientation of the fabric is tilted most, that of the 
deformation gradient is tilted least and thus, the material 
cannot be described by a simple elastic model involving 
only two Lame constants (or bulk modulus and Poisson's 
ratio) as the only parameters. Alternatives are to cut the 
system into pieces with different material properties and 
thus introducing discontinuities jl| or to use the rank four 
stiffness tensor for anisotropic materials [3fJ. The devia- 
toric parts of stress and deformation gradient are seem- 
ingly interconnected via the fabric tensor. 

To conclude, we proposed a consistent averaging for- 
malism to obtain a mean quantity Q in average over arbi- 
trary volumes V. Within this framework, we used for Q c 
the fabric F c = n c <g> n c , the stress er c = (1/V p )f c ® l pc , 
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and the deformation gradient e c = (<5 C / a p )n c <E> n c - A. Fu- 
ture work will involve the extension of the present analysis 
to measure also the fluctuations of the above quantities 
and, in addition, other interesting quantities like, e.g. , 
the plastic strain and non-symmetric parts of the stress 
tensor due to the effective moments acting on single par- 
ticles. More systematic parameter studies are currently in 
progress. 
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